#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
This program draws contour map of the entropy
''' 
from __future__ import division

import common
import os
import numpy
import math
import scipy
import matplotlib.pyplot
import matplotlib.mlab as mlab
import sys

#Returns Entropy for each point that we have temperature
def Entropy(T, E):
    '''
    returns entropy
    '''
    beta = [1 / tt for tt in T]
    S = []
    for i in xrange(len(E)):
        S += [ math.log(2) + E[i] / T[i] + scipy.trapz(E[i:], beta[i:]) ]
        
    minS = min(S)
    return [ts - minS for ts in S]    

#finds two closest points to our T
def closest(T, tList):
    ttList = tList[:]
    ttList.sort()
    for i in xrange(len(ttList) - 1):
        if ttList[i] <= T <= ttList[i+1]:
            return i, i + 1                
    return False
        
#finds weighted average
def average(tx1, tx2, T, tList, S):
     Tmin = tList[tx1]
     Tmax = tList[tx2]
     a = T - Tmin
     l = Tmax - Tmin
     return ( (l - a) * S[tx1] + a * S[tx2] ) / l
            

#matplotlib.pyplot.plot(T, S4, "*", label = "$S = S(\infty) + E(T) / T - \int_{(1/T)}^{0} E(1/T) d (1/T)$")
    
#ttList = os.listdir( os.path.join(path, 'B_' + str(bList[0]), 'L_27', 'Energy') )

maxT = 5
maxB = 5
path = os.path.join("..", "..", "results", "triang", "J_1.0")

bList = numpy.arange(0, maxB, 0.1)
 
Z = numpy.zeros( shape = (maxT * 10-1, maxB * 10 ) )#50 temperature points, 50 magnetic field points

for tb, B in enumerate(bList):
    tList = []
    Energy = []
    result = []
    
    tbPath = os.listdir(path)[:]
    tbPath.sort()
    
       
    for tB in tbPath:
        if round(100 * eval(tB.replace('B_', ''))) == round (100 * B):
            bPath = os.path.join(path, tB, 'L_27', 'Energy')
            break
    
    tList, Energy, err = common.sorted_list(bPath) 
#    for T in os.listdir( bPath ):
#        E = eval( open(os.path.join(bPath, T)).read() )  
#        
#        result += [ (eval(T.replace('T', '')), E) ]        
#    
#    result.sort()
#    
#    tList = [tx[0] for tx in result]
#    Energy = [tx[1] for tx in result]
    
    S = Entropy( tList, Energy )
    
    #roundedB = numpy.around(B, decimals = 1)
    
    for tT, T in enumerate(numpy.arange(0.1, maxT, 0.1)):                 
        #We should figure out what are closest points
        clos =  closest(T, tList)#tx1 and tx2 - indices of the closest temperatures  in the temperature array
        if clos == False:
            print 'ERROR!!!!'       
        tx1, tx2 = clos  
        Z[(tT, tb)] = average(tx1, tx2, T, tList, S)

tList = numpy.arange(0.1, maxT, 0.1)
#bList = numpy.arange(0, maxB, 0.1)


#tList, bList = numpy.meshgrid(tList, bList)
levels = [0.560, 0,480, 0.460, 0.320, 0.240, 0.160, 0.080, 0.060, 0.040, 0.030, 0.020, 0.010, 0.005, 0.001]

CS = matplotlib.pyplot.contour(bList, tList, Z, levels)


B = [0,    0.1,   0.2,   0.3,  0.4,   0.5, 0.6,  0.7,  0.8,  0.9,  1,    1.1,  1.2,  1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2,    2.1,   2.2,  2.3,  2.4,  2.5,  2.6,   2.7, 2.8,  2.9,  3,    3.1,  3.2,  3.3,  3.4,  3.5,  3.6,  3.7,  3.8, 3.9, 4, 4.1,  4.2,   4.3,  4.4,   4.5, 4.6,  4.7, 4.76]
T = [3.62, 3.618, 3.612, 3.61, 3.606, 3.6, 3.59, 3.58, 3.57, 3.56, 3.55, 3.53, 3.51, 3.49, 3.47, 3.45, 3.43, 3.41, 3.37, 3.34, 3.31, 3.275, 3.24, 3.20, 3.16, 3.11, 3.055, 3.0, 2.96, 2.91, 2.85, 2.79, 2.72, 2.65, 2.57, 2.48, 2.4,  2.32, 2.2, 2.1, 2, 1.9, 1.77,   1.61, 1.4,  1.2, 0.9,   0.4, 0]

matplotlib.pyplot.plot(B, T, "*-")
matplotlib.pyplot.xlabel('$B / J$')
matplotlib.pyplot.ylabel(' T / J')
#matplotlib.pyplot.title("Triangular ferromagnetic Ising model in the transverse magnetic field. Phase diagram.")
matplotlib.pyplot.clabel(CS, fontsize=8)
matplotlib.pyplot.show()